Energy spectrum of graphene multilayers in a parallel magnetic field 
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We study the orbital effect of a strong magnetic field parallel to the layers on the energy spectrum 
of the Bernal-stacked graphene bilayer and multilayers, including graphite. We consider the minimal 
model with the electron tunneling between the nearest sites in the plane and out of the plane. Using 
the semiclassical analytical approximation and exact numerical diagonalization, we find that the 
energy spectrum consists of two domains. In the low- and high-energy domains, the semiclassical 
electron orbits are closed and open, so the spectra are discrete and continuous, correspondingly. 
The discrete energy levels are the analogs of the Landau levels for the parallel magnetic field. They 
can be detected experimentally using electron tunneling and optical spectroscopy. In both domains, 
the electron wave functions are localized on a finite number of graphene layers, so the results can 
be applied to graphene multilayers of a finite thickness. 

PACS numbers: 81.05.uf 81.05. ue 73.22.Pr 71.70.Di 
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I. INTRODUCTION 



Graphene monolayers have attracted much attention 
recently because of the unusual Dirac spectrum of elec- 
trons [l|. A remarkable manifestation of the Dirac dis- 
persion is the unusual spectrum of the Landau levels in a 
perpendicular magnetic field, resulting in the anomalous 
quantum Hall effect (QHE) 0-@- These results stim- 
ulated further investigations of the QHE in the deriva- 
tives of graphene. The unusual Landau levels and the 
QHE were obtained for a graphene bilayer in Ref. pj . Al- 
though the Landau levels in graphite were investigated a 
long time ago [8l4l0J , recent studies |llT-[l6| of graphene 
multilayers with a moderate number of layers found in- 
teresting features in the Landau spectrum. Namely, the 
spectrum consists of the two families of levels, whose en- 
ergies scale as B and y/~B, thus indicating the presence of 
both massive and massless Dirac fermions in the system 
[111 ] . The Landau levels for different stacking orders of 
graphene multilayers were studied in Ref. [17| . 

On the other hand, much less attention was paid 
to the orbital effect of a magnetic field parallel to the 
graphene layers. The Shubnikov-de Haas oscillations 
were extensively studied in graphite in a tilted magnetic 
field [n|, but they tend to disappear when the field is 
parallel to the layers. Ref. [18J studied the influence 
of a parallel magnetic field on the putative ferromag- 
netic, superconducting, and metal-insulator transitions 
in graphite. In Ref. [19], the angular magnetoresistance 
oscillations (AMRO) were observed in the stage 2 inter- 
calated graphite (in addition to the Shubnikov-de Haas 
oscillations) for magnetic fields close to the parallel ori- 
entation. AMRO were first discovered in layered organic 
conductors |2Cl [21| and subsequently observed in many 
other layered materials: see, e.g., Ref. [221 and references 
therein. Motivated by the experiment [19f . a theoretical 
study of AMRO in graphene multilayers was the original 
goal of the present paper, with the focus on the peculiari- 
ties due to the presence of two sublattice, Dirac spectrum, 



etc. However, in the standard theory of AMRO [22J, the 

interlayer tunneling amplitude is treated as a small per- 
turbation. It is a reasonable approximation for the inter- 
calated graphite [2J], but not for the pristine graphite, 
where it is generally accepted that the interlayer tunnel- 
ing amplitude is quite large. A non-perturbative treat- 
ment of the interlayer tunneling in a tilted magnetic field 
for graphene multilayers is a very complicated problem. 
So, we decided to focus first on the simpler case of the 
parallel magnetic field. 

Additional motivation for this work comes from the re- 



cent experiment 2J], where the current- voltage I-V re- 
lation was studied for the current perpendicular to the 
layers in a mesoscopic graphite mesa consisting of about 
20-30 graphene layers. The experimental technique is 
similar to the previous work on the cuprate supercon- 
ductors 25j and the charge-density- wave materials |26| . 
When a strong parallel magnetic field up to 55 T is ap- 
plied to the graphite mesa, the dl/dV curve develops a 
peak at a non-zero, magnetic-field-dependent voltage V 
of the order 80 mV. The appearance of the peak may 
indicate formation of the Landau levels in the parallel 
magnetic field, but detailed interpretation of the experi- 
mental results is currently unclear. 

In this paper, we calculate the electron spectrum of 
two or many coupled graphene layers in a strong paral- 
lel magnetic field. To simplify the problem, we consider 
only the minimal model with the electron tunneling am- 
plitudes between the nearest sites in the plane (70) and 
out of the plane (71). The effect of the higher-order tun- 
neling amplitudes [27| is briefly discussed in Appendix [A] 
Our results should be valid for the energies greater than 
the energies of the neglected higher-order tunneling am- 
plitudes and can be verified by tunneling or optical spec- 
troscopy. We focus only on the orbital effect of the mag- 
netic field and disregard possible spin effects [28|. We 
find some mathematical similarities between the electron 
spectrum of graphene multilayers in a parallel magnetic 
field and that of quasi-one-dimensional [29, 30] and quasi- 




FIG. 1: (Color online) A pair of Bernal-stacked graphene lay- 
ers in the parallel magnetic field B applied along the y direc- 
tion. 70 and 71 are the electron tunneling amplitudes. 



two-dimensional [3l| organic conductors |33 |. 

We start with the analysis of a graphene bilayer in 
a parallel magnetic field (Sec. [TTJ) and then proceed to 
the infinite number of layers (Sec. IIII[) . We investi- 
gate both the quasiclassical electron orbits in momentum 
space (Sec. IIIIB Ij) and the exact equation for the en- 
ergy eigenfunctions, which reduces to the Mathieu equa- 
tion (Sec. IIIIB 2]) . We employ both the analytical WKB 
method and exact numerical diagonalization to find the 
energy eigenvalues and eigenfunctions. We identify the 
low-energy domain characterized by closed orbits and dis- 
crete spectrum (Sec. IIIIB 3[) and the high-energy domain 
with open orbits and continuous spectrum (Sec. lIIIB4[) . 
The case of a finite number of layers is analyzed at the 
end of Sec. IIIIB 31 The effect of the tunneling ampli- 
tude 73 responsible for trigonal warping is discussed in 
Appendix [A] 



II. GRAPHENE BILAYER 
A. Model 

First, we consider a graphene bilayer and then gener- 
alize the problem to many layers. The crystal lattice of 
the Bernal-stacked graphene bilayer is shown in Fig. [T] 
The distance between the nearest atoms in graphene 
is a — 1.4 A, and the distance between the layers is 
d = 3.3 A. We restrict our analysis to the minimal tight- 
binding model |16J with the intra- and inter-layer tunnel- 
ing amplitudes 70 = 3.f6 eV and 71 = 0.38 eV. 

There are two sublattices on each graphene layer. 
Thus, the electron wave function is the vector 



v&^V^Vf,^,^), 



(1) 



where the subscripts 1 and 2 enumerate the layers, and 
the superscripts A and B denote sublattices on each 
layer. Sublattices can be selected in various ways. It 
is convenient for us to assign the atoms connected by the 
interlayer tunneling 71 in the Bernal stack to sublattice 
A and other atoms to sublattice B, as shown in Fig.Q] 



E„l 




(a) 




-£ 



e a 




FIG. 2: (Color online) (a) The electron spectrum Q of a 
graphene bilayer in zero magnetic field, (b) The spectrum 
(|12[) in a nonzero parallel magnetic field. The magnetic field 
splits the parabolic spectrum into the two Dirac cones. An 
exaggerated value q — 5 of the magnetic field parameter was 
utilized here. 



In the vicinity of the K point in the Brillouin zone, the 
electron Hamiltonian has the form 



H 



v F (p ■ a) 
liI A 



liI A 
vf(p ■ <f* 



(2) 



Hamiltonian © acts on the vector ((TJ) . Correspondingly, 
cr = (cr x ,(Ty) are the Pauli matrices acting in the sub- 
lattice space; p = p x x + p y y is the in-plane momentum 
measured from the K point; vp = (3/2Ti)"foa « 1 8 cm/s 
is the electron velocity in graphene. The terms vp(p ■ <x) 
and vf{p • <X*) describe the in-plane Hamiltonians of the 
graphene layers. Our choice of the A and B sublattices 
results in the diagonal elements having both <x and cr* 
(complex-conjugated) terms. The term 71/ represents 
the interlayer tunneling, where the matrix 



1 



(I- 



l 




(3) 



connects sublattices A of the adjacent graphene layers. 
Hamiltonian @ has four eigenvalues 



e(p) = ± 



7i 



i \'^ + 4p 2 . 



(4) 



The well-known spectrum ^ is shown in Fig. HJa) . The 
spectrum consists of the four bands with the parabolic 
dispersion for small p. 



B. Parallel magnetic field 



where 



Now let us introduce the in-plane magnetic field B 
yB applied along the y axis. We choose the gauge A 
xBz and use the Peierls substitution 



p — > p+ -A. 

c 



(5) 



Here, we took into account the negative sign of the elec- 
tron charge, so e corresponds to its absolute value. If 
the layer number is denoted by j, the in- plane electron 
momentum on the j'-th layer changes to 



Pj =p + jApx, 



where Ap is 



Ap = -Bd, 
c 



— t = 5 x 10 3 cm _1 T _1 . 
TiB 



(6) 



(7) 



The momentum change Ap has the following physical 
meaning. When an electron tunnels between the layers, 
the Lorentz force F = — -\v x B] changes the in-plane 
momentum by 



Ap x 



F x dt — —B,. 



v z dt = -Bd. 
c 



(8) 



The change in the in-plane momentum results in the rel- 
ative shift of the Dirac points on the different layers in 
the momentum space. 

To simplify equations in the rest of the paper, it is 
convenient to switch to the dimensionless variables 



vpp vpAp e 
>P, > q, > £■ 

7i 7i 7i 



(9) 



Here the parameter q is the dimensionless ratio of the 
"magnetic shift" vpAp and the interlayer tunneling am- 
plitude 71 



vpAp 

7i 



-^ (-Bd) = 0.88 x KT 3 5[T]. (10) 
7i \c ) 



The parameter q describes the orbital effect of the mag- 
netic field in our model and will be frequently referred 
to as the magnetic field for shortness. It is worth noting 
that even for a strong magnetic field this parameter is 
small q < 1, e.g., q = 0.044 for B = 50 T. 

Applying the Peierls substitution ([5]) to Hamiltonian 
@ and switching to the dimensionless variables ©, we 
obtain 



H 



((p-q)-a) 
I A 



I A 

(p-a* 



(11) 



Hamiltonian (|11[) has the following spectrum: 



e(p) = ±-= VP 2 + (P - q) 2 + 1 ± W, (12) 



W 



[{p-q) 2 ~P 2 } 2 + 2p 2 + 2(p-q) 2 + l. (13) 



In contrast to the parabolic dispersion Q, the spectrum 
(TT2l has two Dirac points separated by the magnetic shift 
q with a saddle point in between, as shown in Fig. [2Jb). 
Expanding Eq. (J12I) around one of the Dirac points 



e(p) w ± 



1 



^ 



■.p. 



(14) 



we find that the slope of the Dirac cones is controlled by 
the magnetic field and is greatly reduced for q <C 1 . 

A dispersion similar to Eq. (fi"2l) was found for the 
twisted graphene layers in Ref. [33], where the relative 
displacement of the two Dirac cones in the momentum 
space results from the spatial rotation of the layers 



Ap Iot /h = K A(f> = 6 x 10 6 cm" 



(15) 



Here A<j> = 2° is the twist angle, and Kq = 47r/3\/3a is 
the distance between the Y and K points in the reciprocal 
space. The saddle point between the two Dirac points re- 
sults in the Van Hove singularity in the density of states, 
which was observed experimentally in electron tunneling 
in Ref. [HJ. Comparing Eq. © with Eq. ([15]). we ob- 
serve that the magnetic field effect is much weaker than 
the effect of twisting. Even for B = 50 T, the magnetic 
shift is Ap/h — 2.5 x 10 5 cm -1 is much smaller than the 
rotational shift Ap rot /h. 



III. GRAPHITE 



A. Model 

Now we proceed to the discussion of graphene multi- 
layers. First we solve the problem for an infinite number 
of layers, i.e., for graphite, and then briefly mention the 
effect of a finite number of layers. 

By analogy with the bilayer Hamiltonian ([2]), the 
Hamiltonian of graphite without magnetic field reads in 
the adopted units (O 



/ 



H = 



\ 



(p-a) 


I A 




I A 


(p-er*) 


I A 




I A 


(p-a) 



V 

It acts on the vector 

* = (••• $j-\i>ji>i+x---), 



(16) 



1>i = (4til>f), (17) 



where the subscript j denotes the layer number, and the 
superscripts A and B denote sublattices. 
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FIG. 3: Spectrum ei(p, fc) ((T9J) of Hamiltonian (JTHJ for p y = 0. 
Each curve corresponds to a given value of the out-of-plane 
momentum k indicated on the right. The axes are given in 
the adopted dimensionless units ©. 



Using the momentum representation in the z direction 
and introducing the corresponding momentum k (in addi- 
tion to the in-plane momentum p) , we transform Hamil- 
tonian (|16[) into a 4 x 4 matrix similar to the graphene 
bilayer Hamiltonian ([5} 



H 



(p ■ a) 2I A cos k 
2I A cosk (p-o-*) 



(18) 



Hamiltonian (|18[) has the following spectrum 



ei,2(p,fc) 



±cosfc + \/cos 2 k + p 2 , (19) 

(20) 



£3,4(p, fc) = ± COS fc — v^COS 2 fc 



P~ 



The subscripts (1,2) and (3,4) denote positive and nega- 
tive energies, whereas the subscripts (1,3) and (2,4) cor- 
respond to the terms ± cos fc. Because the spectrum has 
the electron-hole symmetry, we consider only the positive 
energies £1,2- The branches 1 and 2 are equivalent, in the 
sense that £i(p, fc + 7r) = Ei(j>-, k). Thus, it is sufficient to 
consider only one branch £1, which is plotted in Fig. |3] 
Since the off-diagonal elements of Hamiltonian (ITBl van- 
ish for fc = 7r/2, the dispersion has the Dirac-type form 
£i(p, 7r/2) = p for fc = 7r/2, as shown in Fig. [3) 



B. Parallel magnetic field 

1. Semiclassical analysis 

In the presence of a magnetic field, electrons move 
along the isoenergetic surfaces in the momentum space. 
For the field B = yB along the y direction, the elec- 
tron orbits lie on the intersections of the isoenergetic 
surfaces of the dispersion (fTI))) and the planes parallel 
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FIG. 4: Semiclassical electron orbits in the momentum space 
for the in-plane magnetic field along the y axis. Only the or- 
bits with p y — are shown. They are obtained by intersecting 
the {px,k) plane with the isoenergetic surfaces £i(p, fc) = £ 
for (a) e = 0.1, (b) e = 2, (c) e = 2.2. The orbits are (a) 
closed for |e| < 2 and (c) open for |e| > 2. 



to the (p x ,k) plane. The cross-sections of the isoener- 
getic surfaces £i(p, fc) = £ = const with the (p x , fc) plane 
at p y — are shown in Fig. S) The arrows indicate the 
direction of electron motion. 

Topology of the electron orbits changes with the in- 
crease of the energy e. The isoenergetic surfaces for the 
dispersion (TI"9"|) are closed for < e < 2, so the orbits are 
closed too, see Fig. |4ta). Thus, based on the Onsager 
quantization rule [351 ] . the spectrum is discrete for this 
energy interval. However, at the critical energy e = 2, 
the isoenergetic surfaces reconnect, as shown in Fig.@Jb), 
and become open for e > 2, resulting in the open orbits 
shown in Fig. [He). Open semiclassical orbits lead to a 
continuous energy spectrum. 

Fig. 0] shows only the electron orbits for p y — and 
£ > 0. We can find the topology of the electron orbits and 
the character of the spectrum for an arbitrary p y , which 
is a good quantum number for the magnetic field along 
the y direction. The orbits are open, so the spectrum is 
continuous in p x for 



2|£|>p 2 . 



(21) 



The orbits are closed, so the spectrum is discrete and 
degenerate in p x for 



£ 2 -2| £ | < Py <e 2 + 2\e\ 
There are no orbits and no states for 



vl>e 2 



2|£| 



(22) 



(23) 



The domains of the continuous and discrete spectra, de- 
fined by the inequalities (|2"Tj) . (|22p. and (l2"3")l . are shown 
in Fig. [5] in the {p y ,e) plane. 
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FIG. 5: Domains of the continuous and discrete spectra in the 
(p y ,e) plane. The dashed-dotted and dotted curves represent 
solutions of the equations e — 2\e\ — p y and e + 2|e| = 
Py. The spectrum is continuous in the region A defined by 
Eq. (|2ip . discrete in the region B defined by Eq. (I22p , and 
there are no states in the region C defined by Eq. 



2. Mathieu equation 

Now we present a more formal and exact analysis of the 
electron spectrum in a parallel magnetic field. Applying 
the Peierls substitution © in the dimensionless units 



Pj=p + jq, q = xq 

to Hamiltonian (|16[) , we obtain 



(24) 



H = 



(Pj-i ■ t) I A 
I A (Pj ■ <r* 

I A 



I A 
(Pj+i ■ <r) 



V 



■•/ 



(25) 

The eigenvalue problem H^ = £* for H J25J and * (fT7|) 
reads in components 



J A (^-i + V^+i) + 



.(*) 



(p + jq)-e V^=0. (26) 



Here, erM denotes <r for even j and u* for odd j. The 
matrix equation (|26|) represents a set of two equations. 
One of them relates ip? and tp A on the same layer and 
has the simple form 



fl = p x ±ip y +j q ^ 



(27) 



where the signs ± correspond to even and odd j. Using 
Eq. (|27p . we algebraically eliminate ipf components in 
Eq. (|26p and reduce it to the simpler equation 



tf+i+tf-i 



e — 



(p + jq) 



ti 



(28) 



which has the same form for even and odd j. From now 
on we drop the superscripts A. In the Fourier represen- 
tation 



Eq. ((551) becomes 

d_ _ .Px 
dk a 



il)(k)e lk3 dk, 



i>(k)-v(k)i>(k) = o, 



where 



V(k) 



2s 



cosfc 



(29) 



(30) 



(31) 



Here, ip(k) is a 27r-periodic, twice-differentiable function 
ijj(k) = ijj(k + 2ir). To further simplify Eq. ([30]). we in- 
troduce the function <j){k) 



^(fc)=e 2fe( ^/«V(fc), 



(32) 



which eliminates the term ip x /q from Eq. (|30[) and re- 
duces it to the angular Mathieu equation for </>(fc) 



d 2 (/>(k) 
dk 2 



V(k)cj)(k) = 0. 



(33) 



Eq. p3p is equivalent to the Schrodinger equation for 
a particle moving in the ID potential V(k) (|3ip . The 
variables e and p y are the parameters that control V(k). 
Since V(k) is periodic in k, the Bloch theorem can be 
applied, so the solutions of Eq. ((33]) have the form 



Mk) 



*u K {k). 



(34) 



Here, k is the quasimomentum in the space reciprocal 
to the k space, and u K (k) is a 27r-periodic function in k. 
From Eqs. f|32|) and (|34p and the periodicity requirement 
for i/j(k), we select the solutions of Eq. f|33|) with k — 
—p x /q. Since the solutions of Eq. (I33P are periodic in 
the quasimomentum cj) K +i(k) = <p K (k), the parameter e 
in Eq. (|30p must be periodic in p x : e(p x ) — e(p x + q). 
Therefore, the magnetic field effectively introduces the 
magnetic Brillouin zone in p x with the period q. 

We have reduced the original eigenvalue problem ((25)) 
to the convenient differential equation ((SB"]) . Different 
regimes for its solutions are controlled by the parameters 
p y and e. If the criterion (|22|) is satisfied, the ID classical 
motion is bounded by the barriers of V(k), as shown in 
Fig. (6^a) for e = 0.1 and p y = 0. Then the energy spec- 
trum is discrete. On the other hand, if the criterion (|2ip 
is satisfied, the potential is negative V(k) < for any 
k, as shown in Fig. |U(b) for e = 2.1 and p y = 0. Then 
the motion of a particle is unbounded, and the spectrum 
is continuous in p x . The first regime corresponds to the 
closed orbits in Fig. QJa), and the second regime to the 
open orbits in Fig. [He). In the following sections, we use 
the approaches of both Sec. IIIIB II and this section to 
obtain and interpret the results. 
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FIG. 6: The plots of V(k) (O for p y = 0. The cases of 
e = 0.1 and e = 2.1 are shown on the panels (a) and (b). 
The classically permitted region corresponds to V(k) < 0, as 
indicated by the thick horizontal line. Thus, the panels (a) 
and (b) represent bounded and unbounded motion. 



3. Closed orbits 

In this section, we study the electron spectrum in the 
domain of the (p y ,e) plane defined by the criterion ([22)) 
and labeled by the letter B in Fig. [5) In this case, the 
classical motion of a particle in the ID potential (|3"Tj) is 
restricted to the potential wells separated by the barriers, 
as shown in Fig. UJa). The height of the barriers is 



h(p. 



!)• ' 



max[V(k)] 



e(2 



Pi 



<r 



(35) 



The barriers h(p y ,e) 3> 1 are high everywhere, except at 
the boundary of the domain (1221) . Thus, we can neglect 
tunneling and use the WKB quantization rule for a single 
well of V(k) 



-2e cos k + e 2 — p\ dk = 2tt [ n + — ) q, 



arccos a 



where 



2 2 

£ -Py 
2e 



(36) 



(37) 



The integral on the left-hand side of Eq. (|3"6l is propor- 
tional to the area enclosed by the electron orbit in mo- 
mentum space, see Fig.@|a). Thus, Eq. (|36[) is equivalent 
to the Onsager quantization rule in a magnetic field 35] . 
Using the incomplete elliptic function of the second kind 



E{4>,m) 



vT 



m 2 sin a da, 



(38) 
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FIG. 7: (Color online) Low-energy levels e n (px,p y ) for q = 
0.044. Panel (a) shows e n (p x ,0) vs p x for p y — 0, and panel 
(b) shows e„(0,p y ) vs p y for p x = 0. Solid lines represent 
exact numerical diagonalization of Hamiltonian (|25[) . Dashed 
lines represent the WKB analytical approximation (|39|) . All 
quantities are presented in the adopted dimensionless units 



Eq. (|36|) can be written as 



8V2eVTT^E 



7r + 2 arcsin a 



l + o 



2tt ( n + - ] q. 



(39) 



Eqs. ([39| and ([37)) implicitly define s n (p x ,p y ) as a func- 
tion of p y and n for a given magnetic field q, and the 
spectrum is degenerate in p x . 

To check validity of the WKB approximation, we diag- 
onalize of the original Hamiltonian (|25j) numerically and 
compare results with the solutions of Eq. (f3"!)f . Momen- 
tum dependences of e n (p x ,0) and e n (0,p y ) are shown 
in Figs. E{a) and (b) for a few lowest energy levels at 
q = 0.044. The analytical approximation (|39l) agrees 
well with the numerical results for n/ 0. The discrete 
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FIG. 8: (Color online) Energy levels e n (0,0) vs the quantum 
number n at p x = and p y = for q = 0.044. The horizontal 
axis shows the combination q(n + 1/2). The circles represent 
solutions of the WKB equation (|39[) . and the small points 
inside the circles represent numerical data. The quadratic 
approximation ()41[) is shown by the solid line. All quantities 
are presented in the adopted dimensionless units Q. 



energy levels shown in Fig. EJa) are degenerate in p x and 
represent the Landau levels in a parallel magnetic field. 
However, the n = level has a remarkable dispersion in 
p x . Similarly to the spectrum of the graphene bilayer in 
Fig. EJb), the n = level consists of a series of the Dirac 
cones shifted by the vector q. This dispersion cannot be 
obtained from the approximate WKB equation (|39| , be- 
cause of the divergence at e = in the original equations 
(EU and (USD. 

Fig. [7jb) shows that the energy levels e n (0,p y ) have 
a quadratic dispersion in p y , except for the n = level. 
Given the degeneracy in p x , the energy levels £ n (Px,Py) 
form one-dimensional bands in p y , so the density of states 
diverges at the bottom points e n (p x ,0) of the bands. 
These singularities in the density of state can be de- 
tected experimentally by electron tunneling or optical 
spectroscopy. The energies e ra (0, 0) are plotted in Fig. [5] 
in the interval < e < 2 vs the combination q(n + 1/2) 
appearing in Eq. (|39|) . Depending on the magnetic field 
q, a different number n max of the discrete levels fills the 
curve. By setting e = 2 in Eq. (j3"6"j) . we obtain 



1 _ 8\/2 _ 3.6 

2 irq q 



(40) 



For example, for q = 0.044, we have n max = 81 levels, 
which are depicted by circles in Fig. [5] 

For small e and p y — 0, we find from Eq. (|39[) 



0.7 g* I n+- 



(41) 
We observe that the energies (|4"T1) depend quadratically 
on the level number n and the magnetic field q. This 
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FIG. 9: (Color online) The absolute values of the wave func- 
tions \ipf\ an d \ipf\ on the sublattices A and B vs the layer 
number j for q — 0.044, p x = q/2, and p y = 0. The variable 
n denotes the energy level number. 



dependence is different from the usual Landau level de- 
pendence, where the energies are linear in the field and 
in the quantum number. The reason for the unusual de- 
pendence in our case is the following. For small e, the 
semiclassical orbit in Fig. Ufa) shrinks to a thin ellip- 
soid of the length 7r in the k direction and the width 
V2~£ in the p x direction. Thus, the area enclosed by the 
semiclassical orbit is proportional to y/e, so the Onsager 
quantization rule gives quadratic dependence of the en- 
ergy on the magnetic field and the level number n. The 
quadratic approximation (|41 [) is shown by the solid line 
in Fig. |8] and works well in the region e < 0.1. 

Fig. |9] shows the plots of \ipf\ and \ijjf \ vs the layer 
number j for several energy levels n. We observe that 
the magnetic field causes localization of the wave function 
on a finite number of layers. According to Eq. (|27[) . the 
wave functions for the low energy levels e n are localized 
predominantly on the sublattice B. The magnitudes of 
|^| and \tpf\ are shown by different vertical scales in 
panels (a) and (b) of Fig. [9] 

Now let us briefly discuss the spectrum of a finite sys- 
tem with the total number of layers N. Fig. [TU] shows 
£ n (Px, 0) for N = 7, 21, and oo. The degeneracy in p x is 
lifted for a finite number of layers, but, with increasing 
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FIG. 10: (Color online) Energy spectrum e n (px,0) vs p x at 
p y = and q — 0.044 for the system with a finite number of 
layers N. The dashed-dotted, dashed, and solid lines corre- 
spond to N — 7, 21, and oo. All quantities are presented in 
the adopted dimensionless units (|9}. 



N, the spectrum approaches to that of the infinite sys- 
tem with N — oo. Indeed, if the localization length for 
a particular energy level is shorter than the size of the 
system, the energy of the level is the same as for N = oo. 
Thus, the results obtained for N = oo are applicable to 
a finite system with a sufficient large N. 



4- Open orbits 

Now we study the energy spectrum in the domain de- 
fined by Eq. (l2l"j) and labeled by the letter A in Fig. [5] 
It corresponds to the open electron orbits in Fig. Etc). 
In the Mathieu equation (|55)l . the potential V(k) < is 
negative for any k, so the motion is unbounded, as shown 
in Fig.Etb). Then, the WKB solutions of Eq. ([33]) are 



4>{k) 



= e ±iS(k) 



\ S(k) = / VWW\dk, (42) 



where the signs ± correspond to the direction of mo- 
tion. Because of the periodicity requirement for ip(k) 
and Eq. (|3"2"|) . the phase accumulation in Eq. f4"2")l over 
the period 2n must be equal to —2i:p x /q plus an integer 
multiple of 2ir. Thus we obtain the following quantiza- 
tion condition for the open orbits 



gS*(27r) 



-2ecosfc + e 2 



Pydk 



T^(Px + nq). 

Here the integer n is different from the integer n in 
Eq. (l3l)|) and takes both negative and positive values. 
The sign of p x + hq corresponds to the two solutions in 



FIG. 11: (Color online) Energy spectrum around e — 2 for 
q = 0.044 and p y = 0. The solid lines are obtained by numer- 
ical diagonalization of Hamiltonian (|25l) . The dashed lines, 
obtained from Eq. (|39|1 . are labeled by the integer n shown 
on the right. The dashed-dotted lines, obtained from Eq. (|44[) . 
are labeled by the integer h shown at the top. This plot illus- 
trates a transition from discrete to continuous spectrum. All 
quantities are presented in the dimensionless units @. 



Eq. (j42|) . Eq. (|43|) can be represented in terms of the 
elliptic integral (f38|) 



4V2eVl + a E 



2'Vl + fl 



= T^{Px + nq), (44) 



where the parameter a is given by Eq. (|37p. In contrast 
to Eq. (|39l) . Eq. (|44|) contains p x , so the energy levels 
£ n(Px,Py) continuously depend onp^. 

The energy spectra obtained from Eqs. ([3T)]) and ([4^]) 
are compared in Fig. [11] with the results of numerical di- 
agonalization of Hamiltonian (|25[) around the critical en- 
ergy e — 2 for p y = 0. For e < 2, the spectrum consists 
of the energy levels degenerate in p x , which are well de- 
scribed by the analytical approximation (|39p for closed 
orbits. The corresponding level number n is shown on 
the right in Fig. [Tl] At the energy e = 2, the spectrum 
undergoes a transition to the regime of continuous dis- 
persion uip x . For e > 2, the spectrum consists of the two 
families of lines with the opposite slopes. This spectrum 
is well described by the analytical approximation f4"4"]) for 
open orbits. The corresponding number n labels the dis- 
persion lines and takes both positive and negative values 
shown at the top in Fig. [11] Because the left-hand sides 
of Eqs. ((351) and (O differ by the factor of 2 at e = 2 and 
p x = 0, the numbers n and n are connected asnw 2|n|. 
The approximations ([35]) and (j4"4"]) stop working in the 
vicinity of the critical energy e = 2. 

For a high energy e, when the parameter a (|37p is large, 
we can obtain the spectrum explicitly by expanding the 
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FIG. 12: (Color online) The absolute value of the wave func- 
tion \ipj\ vs the layer number j for a state from the do- 
main of continuous spectrum. The parameters of the plot are 
n = -250, En = 10.98, p x = p y = 0, and q = 0.044. The solid 
and dashed lines represent exact numerical diagonalization 
of Hamiltonian ((25} and the approximate analytical formula 
([50)1 . respectively. 



square root in Eq. (|42)l for S(k) in powers of 1/a 

e sin k 



Pi. 



<7\/£ 2 -pl 



S(k) 



Then the quantization condition (|43| gives 
ye 2 -p 2 y = T(Px + nq), 



which can be written as 



el 



(Px + nq) 2 + p 2 



(45) 



(46) 



(47) 



The spectrum Eq. (|4"T)) is the same as for decoupled 
graphene layers with the Peierls substitution ([24]) . 

Substituting the square root expression from Eq. (|46[) 
into Eq. (l4"5l) . we obtain the approximate WKB wave 
functions H2} for e > 



4>{k) = exp 



.Px + nq . 
-I k + i 



d \Px + nq) 2 + p 2 y 



q(j>x + nq) 



sinfc 



(48) 

Then, using Eq. (1521) , we calculate the Fourier transform 
(|29p and find the wave function ipj in the direct space 



ipj = Jh-j 



^(Px + nq) 2 +pl 
q(p x + nq) 



(49) 



Here J m (x) is the Bessel function of the ra-th order of 
the first kind. For p y — 0, Eq. (l4"9l simplifies to 



Ipj = Jf L - 



sign(p x + nq) 



(50) 



The wave function (|50j) is centered at j = h, as shown 
in Fig. [T2J We observe that, even though Eq. (14"T1) coin- 
cides with the spectrum of effectively decoupled graphene 
layers, the corresponding wave function (|50[) is localized 
on a large number of layers proportional to 1/q. Similar 
wave functions are known for the quasi-one-dimensional 
conductors in a magnetic field [29|, [3CJ • 



IV. CONCLUSIONS 

In this work, we have studied the orbital effect of a 
strong magnetic field applied in the y direction parallel 
to the layers of the graphene bilayer and multilayers. For 
the former, the magnetic field splits the parabolic bilayer 
dispersion into the two Dirac cones in the momentum 
space with the spacing proportional to the magnetic field. 
For the latter, we have found two domains in the parame- 
ter space with distinct energy spectra. In the low-energy 
domain, the scmiclassical electron orbits are closed, so 
the spectrum is discrete and degenerate in p x . The en- 
ergy levels depends quadratically on p y , thus forming a 
series of one-dimensional bands in p y . The discrete en- 
ergies of the bottoms of the bands are the analogs of the 
Landau levels but depend quadratically on the energy 
level number n and the magnetic field B. The n = en- 
ergy level around zero energy has unusual properties and 
consists of a series of shifted Dirac cones, similarly to the 
bilayer case. In the high-energy domain, the semiclassical 
electron orbits are open, so the spectrum in continuous in 
p x , thus forming two-dimensional bands in p x and p y . For 
high enough energies, these bands evolve into the Dirac 
cones originating from different layers and shifted in the 
momentum space due to the applied magnetic field. In 
both regimes, the wave functions are localized on a finite 
number of layers. Mathematically, the problem reduces 
to the Mathieu equation. The WKB approximation for 
the semiclassical electron orbits in the momentum space 
in the magnetic field agrees well with exact numerical di- 
agonalization, except for a few special cases, where the 
WKB approach is not applicable. 

The obtained energy spectrum can be verified experi- 
mentally using electron tunneling or optical spectroscopy. 
Our results may help to understand the experimentally 
measured I-V curves for a mesoscopic graphite mesa in 
a strong parallel magnetic field up to 55 T [2J| , although 
detailed interpretation is not clear at this point. 

We studied the minimal model with the two tunnel- 
ing amplitudes between the nearest neighboring sites in 
the plane (70) and out of the plane (71). In general, 
the obtained results should be valid for the energies 
greater than the neglected higher-order tunneling am- 
plitudes [271 ] . which can be taken into account in future 
studies, if necessary. A more detailed discussion of the 
influence of the trigonal warping amplitude 73 on our 
results is given in Appendix [A"l 
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Appendix A: The effect of the higher-order 
tunneling amplitudes 

In this appendix, we examine how our results are af- 
fected by inclusion of the higher-order tunneling am- 
plitudes [23]. We focus on the tunneling amplitude 
73 = 0.29 eV [27j . which connects any given atom B 
with the three nearest atoms B on the adjacent layers 
(see Fig. Q}. This term has the C3 rotational symmetry 
and is responsible for the trigonal warping of the electron 
spectrum. Given that the amplitudes 73 and 71 are of the 
same order, we need to explain why one can disregard 73, 
but use 71 at the same time. 

In the presence of 73, the model Hamiltonian becomes 



H = 




(Al) 



where c.c. means complex conjugated, and we introduced 
P = Px + ip y and the small dimensionless parameter 



73 



73 

To 



0.08. 



(A2) 



Because the 73 term vanishes at the K point (p = 0), 
we expanded this term in Eq. (|A1|) to the first order in 
p, to be consistent with the linearization of the 70 term. 
In contrast, the 71 term does not vanish at the K point. 
Thus, the 73 term in Eq. (|Alj) is smaller than the 71 
term by the small factor 73 for the energies where the 
linearization in p is applicable, even though 71 and 73 
are of the same order. For this reason, the 73 term can 
be generally neglected relative to the 71 term, except for 
the very low energies, as discussed below. 

Hamiltonian (IA1[) is written in the conventional system 
of units. In the system of units (0), adopted in our work, 
Hamiltonian becomes 



H 



p 2cosfc 

p* 2j' 3 pcosk 

2 cos A: p* 

2j' 3 p*cosk p 



(A3) 



Hamiltonian (|A3|) differs from Hamiltonian ([18|) by the 
additional terms proportional to 73. 

Our results for the electron spectrum in a parallel mag- 
netic field rely mainly on the shape and topology of the 



isoenergetic surfaces e(p, k) = e = const in the momen- 
tum space. The isoenergetic surfaces for Hamiltonians 
(|A3|) and (fT8")) are compared in Fig.[T3]for e — 0.2 in panel 
(a) and for e = 3 in panel (b). The green, cylindrically- 
symmetric surfaces correspond to Hamiltonian (|18[) (their 
cross sections are shown in Fig. 2]), whereas the red 
surfaces with the C3 symmetry correspond to Hamilto- 
nian (IA3|) . We observe that 73 produces trigonal warp- 
ing of the isoenergetic surfaces, but the topology of the 
surfaces does not change for |e| > 0.1. Thus, our con- 
clusions about the discrete energy spectrum for |e| < 2 
and continuous for |e| > 2 in the presence of a paral- 
lel magnetic field remain qualitatively valid. However, 
the energy spectrum acquires anisotropy with respect to 
the in-plane rotation of the magnetic field, which can be 
studied experimentally. 

Trigonal warping of the isoenergetic surfaces becomes 
progressively more pronounced at low energies |e| < 0.1. 
At much lower energies e ~ 0.01, the isoenergetic surface 
spits into four separate branches with three new Dirac 
points surrounding the original Dirac point [lOJ 1271 ] . For 
such low energies, the isoenergetic surfaces qualitatively 
differ from the case of 73 =0, and our results are not ap- 
plicable. At the low energies, it is also necessary to take 
into account the other tunneling amplitudes, in particu- 
lar 72 ~ 10 meV, which connects the next-nearest layers 
[271 ] . In the presence of many tunneling amplitudes, the 
problem becomes very complicated. In addition, disorder 
may smear out the spectrum at low energies. 

Thus, we restrict the applicability of our results to 
the relatively high energies |e| > 0.1, which is about 
40 meV in the dimensional units. The energy spectrum 
in this range can be studied by tunneling or optical spec- 
troscopy. Notice that the magnetic-field-dependent peak 
in dl/dV was observed in Ref. 24| at the applied voltage 
of about 80 meV, although exact nature of this peak is 
still unclear. 




FIG. 13: (Color online) Isoenergetic surfaces e(p x ,p y ,k) = 
e — const are shown for e = 0.2 and e = 3 in panels (a) 
and (b), respectively. The green cylindrically-symmetric sur- 
faces correspond to Hamiltonian (|18|) . and the red trigonally- 
warped surfaces to Hamiltonian (|A3|) . 
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